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ABSTRACT 

We present new interferometer observations of the ^^CO (1-0) and ^^CO 
(2-1) line emission of NGC 1068 with a resolution of 0.7". The molecular gas 
in the inner 5" is resolved into a ring with two bright knots east and west of 
the nuclear continuum emission. For the first time in NGC 1068, we can trace 
molecular gas at ~ 0.18" (13 pc) from the nucleus. The high velocities in this 
region imply an enclosed mass of ~10^ Mq. This value is consistent with a 
black hole mass of 1.7 x 10^ M0, as estimated from nuclear H2O maser emission, 
plus a contribution from a compact nuclear stellar cluster. Perpendicular to the 
kinematic major axis optical images of NGC 1068 show a bright, stellar, oval 
structure of eccentricity 0.8 and a deprojcctcd length of 17 kpc. Analysis of the 
rotation curve shows the CO spiral arms are at the inner Lindblad resonance of 
this bar-like structure. Inside the molecular spiral arms, 10" from the nucleus, 
the CO kinematic axis changes direction probably in response to the to the 
2.5 kpc (deprojected) long stellar bar seen in the near infrared (NIR). The low 
velocity dispersion indicates the molecular gas is in a disk with a thickness of 
10 pc in the nuclear region and 100 pc in the spiral arms. 

We constructed kinematic models for the molecular gas using elhptical orbits 
caused by a ~ 1" (72 pc) nuclear bar and using tilted rings resulting in a warp. 
We find that the gas motions are consistent with either the warp or the bar 
models. However, because there is no evidence for a ~ 1" nuclear bar in NIR 
images, we favor the warp model. A warped CO disk can also explain the 
obscuration of the AGN, the extinction of light from the nuclear stellar cluster, 
and the observed NIR and mid-IR polarization. The model predicts the warped 
CO disk should become edge-on at a radius of 70 pc, thereby creating a cavity 
for the ionization cone. 



Subject headings: galaxies: ISM - galaxies: nuclei of - radio lines: ISM - 
galaxies: individual (NGC 1068) 
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1. INTRODUCTION 

Knowledge of the distribution and kinematics of the circumnuclear molecular gas in 
Seyfert galaxies is essential for understanding the fueling of the central region and the 
role of gas and dust in obscuring the active galactic nucleus (AGN). Recent improvements 
in millimeter interferometry now allow us to observe molecular lines with sub-arcsecond 
resolution, high spectral resolution, and high sensitivity. In nearby Seyfert galaxies, with 
the best resolution currently available, we can now search for molecular gas down to radii 
of a few tens of parsecs. Models of the motions of this molecular gas can then help us to 
understand the mechanism that transports gas into the center. 

An overview of the modeling results so far is given in Schinnerer, Eckart & Tacconi 
(1999a; hereafter SET99a), and a detailed model of the nuclear molecular gas in the 
Seyfert 1 galaxy NGC 3227 is given by Schinnerer, Eckart & Tacconi (1999b; hereafter 
SET99b). In the inner 80 pc of NGC 3227, for example, the kinematics of the molecular gas 
are better explained by a warping of the thin gas disk than by motion in the potential of a 
nuclear stellar bar. In that galaxy, molecular gas is detected down to a radius of ~ 10 pc. 

In this paper, we model the complex kinematics of the central few hundred parsecs of 
the Seyfert 2 galaxy NGC 1068. This galaxy is at a distance of 14.4 Mpc [Hq = 75kms~^; 
see Bland-Hawthorn et al. 1997), giving a scale of 72 pc arcsec"^ (Table which is very 
favorable for detailed studies of the circumnuclear gas. NGC 1068 was originally classed as a 
Seyfert 2 galaxy based on its direct-path, i.e. unpolarized emission lines, which have narrow 
widths (Khachikian & Weedman 1974). This narrow-line region (NLR) is extended to the 
northeast of the nucleus and coincides with the radio jet. Antonucci & Miller (1985) later 
showed that the scattered-path, polarized emission lines have widths typical of broad-line 
regions (BLR), indicating that NGC 1068 has a hidden Seyfert 1 nucleus. Around this 
nucleus, there is a central nuclear star cluster with a FWHM diameter of 0.6" (43 pc) that 
was mapped by Thatte et al. (1997). On a much larger scale, there is a stellar bar 2.3 kpc 
long, detected in the NIR (Scoville et al. 1988; Thronson et al. 1989). This bar is in turn 
surrounded by a circumnuclear starburst ring that, along with other star forming regions, 
contributes half the bolometric luminosity of the galaxy (Telesco et al. 1984; Davies et al. 
1998). This starburst ring is also detected in lines of molecular gas (Planesas et al. 1991; 
Jackson et al. 1993; Tacconi et al. 1994; Heifer & Blitz 1995; Tacconi et al. 1997) and in 
the radio continuum (e.g. Gallimore et al. 1996). 

The plan of this paper is as follows. Section 2 describes the observations and data 
reduction. Section 3 covers the distribution and kinematics of the CO emission, the 
millimeter continuum, and the location of the nucleus. In section 4 we derive the rotation 
curve of the molecular gas and discuss its relation to the stellar bar and the location of 
dynamical resonances. We then derive the molecular gas mass and the total dynamical mass 
(section 5) and derive the thickness of the molecular gas disk from the velocity dispersion 
(section 6). In section 7, we fit the kinematic data with bar and warp models, and discuss 
the model results in section 8. Section 9 summarizes our conclusions. 
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2. OBSERVATIONS AND DATA REDUCTION 

The ^^CO lines were observed in February/March 1997 with the IRAM interferometer 
on Plateau de Bure, France, with 5 antennas in the long-baseline A, Bl, and B2 
configurations and in November 1995 to March 1996 with 4 antennas in the shorter-baseline 
CI, C2, and Dl configurations. The ^^CO (1-0) observations had baselines from 24 to 
408 m that gave a beam of 2.2" x 1.2" at p. a. 26° after uniform weighting. The ^^CO (2-1) 
observations had only the A, Bl, and B2 baselines, ranging from 40 to 408 m, that gave 
a beam of 1.0" x 0.5" (p. a. 26°) after uniform weighting. The bandpass was calibrated 
on 3C454.3 and 3C345. The source 0235+164 was observed every lOmin for phase and 
amplitude calibration. The data were CLEANed with 700 iterations and convolved to 
beams of 1.4" x 1.4" at ^^CO (1-0) and 0.7" x 0.7" at ^^CO (2-1) , with velocity channels 
of lOkms"^ for both lines. We also made a 115 GHz continuum map from the line-free 
channels from —310 to — 190kms~^ and +200 to +320kms~^. To estimate the line fiux, we 
subtracted the mean continuum from the channels with (line+continuum). 

Our ^^CO (1-0) map agrees well with CO maps made with other interferometers, and 
our integrated ^^CO (1-0) fiux agrees with that in the interferometer spectrum of Heifer 
& Blitz (1995) to within 15%. Heifer & Blitz (1995) also combined their BIMA array and 
NRAO single dish data. Their BIMA/NRAO map has 25 to 30% more fiux than their 
BIMA map alone, so some of the large-scale fiux may be resolved out by the BIMA and 
IRAM arrays. To check our calibration of the ^^CO (2-1) line emission, we made a ^^CO 
(2-1) map with the same beam as our ^^CO (1-0) map. At a resolution of 1.4", the qo|^Io| 
ratio for the nuclear ring is 1.0 (±20%), in agreement with the ratio of 1.1 found in this 
region with the IRAM 30 m telescope (Braine & Combes 1992). 



3. MOLECULAR GAS DISTRIBUTION AND PROPERTIES 

3.1. The distribution of the line emission 

Our new observations allow us to study the millimeter molecular line and continuum 
emission at an unprecedented angular resolution. The ^^CO (2-1) data, in particular, give 
very detailed maps of the molecular gas in the inner 200 pc of NGC 1068. At 100 pc from 
the center, the ^^CO (2-1) emission is in a ring (hereafter the nuclear ring; Fig. |l]). In 
this nuclear ring, there are two emission knots, one at 1" east of the nucleus and another 
one 1.5" west and 1" south of the center. At ±0.5" north and south of the center, there 
are emission ridges connecting the two knots at 25% of the peak line intensity. This same 
overall nuclear structure also exists in the CS(2-1) line (Tacconi et al. 1997) and the 
HCN(l-O) line (Tacconi et al. 1994). The go|^Io| ratio is between 0.8 and 1.0 over the 
entire nuclear ring, except for a part of the western emission knot, which has a ratio of 2.0. 
This is also evident in the nuclear spectra. The ^^CO (2-1) profile has three peaks while 
the ^^CO (1-0) profile has only two. The missing peak in the ^^CO (1-0) line is offset by 
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+75kms ^ relative to the systemic velocity and is associated with the western knot. This 
suggests the CO lines are optically thin in parts of the western knot. 

As shown in earlier maps, most of the ^^CO (1-0) emission is in two spiral arms with a 
diameter of ~ 40" (Fig. |l|; Planesas et al. 1991; Kaneko et al. 1992; Heifer & Blitz 1995). 
On our new map, the ^^CO (1-0) peak ~ 1" east of the center is 1.5 times stronger than the 
spiral-arm emission. In our map, the northern spiral arm is less luminous than the southern 
one, in contrast to the earlier measurements of Heifer & Blitz (1995; see below). 

CO emission is also observed along the NIR stellar bar found by Scoville et al. (1988) 
at about half of the nuclear peak intensity. With our 1.4" beam, we mainly find emission in 
the northern half of the bar. In the lower resolution (~ 3.4") ^^CO (1-0) data (Tacconi et 
al. 1997) emission is seen on both sides of the bar, but stronger north of the nucleus. Both 
sides of the bar were also detected in lower resolution ^^CO (1-0) maps by Heifer & Blitz 
(1995) and in the ^^00(1-0) map by Tacconi et al. (1997). The tips of the NIR stellar bar 
almost touch the spiral arms. The spiral arms are not detected in the ^^CO (2-1) map, 
because their diameter is larger than the 22" primary beam of the IRAM antennas, and a 
mosaic map would be needed to recover the spiral arms in ^^CO (2-1) . Furthermore, the 
spiral arms are weaker than the nucleus and they may be partly resolved since we only used 
long baseline configurations for the ^^CO (2-1) mapping. 

Our ^^CO (1-0) intensity map agrees well with the OVRO map (Baker & Scoville 1998; 
private communication). Differences with the earlier data of Heifer & Blitz (1995) may be 
due to uncertainties in the calibration and deconvolution of the BIMA/NRAO data and the 
possibility that 30% of the extended line flux is resolved out in our interferometer map. 



3.2. Radio continuum emission 

Our 115 GHz continuum map shows two emission knots, one coinciding with the 
nucleus, and the other 2.5" east and 4.5" north of the nucleus (Fig. |^). The nuclear knot 
is slightly resolved by our 1.4" beam. A Gaussian fit gives 1.5" x 1.7" at p. a. 156°. At 
115 GHz, the nuclear knot has a flux density of (36±5) mJy, while the northern knot has 
(15±5)mJy. The northern knot coincides with the peak of the 5 GHz radio jet mapped 
by Wilson & Ulvestad (1983), as noted previously by Heifer & Blitz (1995), who observed 
the same two continuum knots in their BIMA map at 85 GHz. This structure at 85 GHz is 
also observed with the IRAM array (Tacconi 1999; private communication). At 85 GHz, 
both knots have the same flux density, 40 mJy. At 230 GHz, no continuum emission was 
detected, to a limit of 6mJy at each knot. 
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4. GLOBAL KINEMATICS OF THE MOLECULAR GAS 

In this section, we describe the velocity field of the molecular gas, estimate the enclosed 
nuclear mass, and derive a rotation curve. We then combine this rotation curve with optical 
and IR evidence for bars in NGC 1068 to deduce the locations of dynamical resonances. 



4.1. The large-scale CO velocity field 

A striking feature of the ^^CO (1-0) velocity field is the change of position angle of the 
kinematic major axis within the spiral arms (Fig.|^ a and c). At radii > 15" the velocity 
field of the molecular gas joins the outer HI velocity field at p. a. 280° and inclination 40° 
(Brinks et al. 1997). At r < 10", the kinematic major axis changes its position angle from 
278° to 305°. The change in position angle at r ~ 10" also occurs in the velocity fields of 
the Ha (Dehnen et al. 1997) and the stars (Garcia-Lorenzo et al. 1997). This effect also 
appears in the HCN(l-O) line (Jackson et al. 1993; Tacconi et al. 1994) and was analyzed 
by Heifer & Blitz (1995) in ^^CO (1-0) . 



4.2. The nuclear kinematics from the CO position-velocity diagrams 

In the inner 500 pc, the ^^CO (2-1) velocity field is perturbed by the two bright knots 
~ 1" east and west of the nucleus that do not follow the overall circular rotation. To more 
easily separate symmetric and asymmetric components and see which bits of the molecular 
gas have ordered motions, we used position-velocity diagrams {pv diagrams) along different 
angles. We made pv diagrams every 20° starting at 10° (Fig. The diagram along p. a. 
110° is remarkably symmetric and has most of the flux, because it crosses both emission 
knots. All diagrams close to the kinematic minor axis are compact, because there is only 
weak emission north and south of the center. 

In ^^CO (2-1) , the western knot (knot W) is slightly farther from the nucleus (~ 1.5") 
than the eastern knot. Knot W appears in pv diagrams at p. a. 50° and 70°, with a velocity 
dispersion of 60kms~^. Knot E, with a velocity dispersion of 70kms~^, appears in pv 
diagrams at p. a. 150°, 170°, 0°, and 10° as a distinct component from —230 to —30 kms~^ 
at ~ 1.5" southeast of the center. To estimate these components' flux, we calculated the 
moments only for velocities —230 to — 30kms~^. The apparent dispersion of knot W is a 
mixture of the 40kms~^ velocity spread of the nuclear disk with the 30kms~^ spread of the 
knot itself. For the eastern knot, it is hard to see the knot component, since it is blended 
with the larger-scale rotating disk. These two knots (Table |2D contain a few 10^ Mq of 
molecular gas, which suggests they are giant molecular cloud complexes. There is a similar 
mixture in the circumnuclear molecular gas in NGC 3227, where we also found several 
components deviating from the overall circular rotation (SET99b). 
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4.3. Hints of a nuclear enclosed mass 

All position-velocity diagrams show significant high velocity emission at r < 0.2", 
similar to that in the nucleus of NGC 3227 (SET99b). If the high velocities trace rotation 
in the immediate vicinity of the nucleus, we can estimate the enclosed nuclear mass. The 
minimum radius at which we detect an emission peak of high-velocity nuclear molecular 
gas is 13.5 pc, and its mean velocity relative to systemic is Ai; ^ 188kms~^, uncorrected 
for inclination. If this is rotation, it implies a mass > 1.5 x 10^ Mq within the inner 27pc 
diameter, comparable with the mass of 4.5 X 10^ Mq within the inner 70 pc of our Galaxy 
(Genzel et al. 1994). This is ten times the black hole mass of 1.5 x 10^ estimated from 
H2O masers in the inner 2pc of NGC 1068 (Greenhill & Gwinn 1997). The mass in the 
inner 27pc is thus mostly due to a central, compact part of the 0.6" (43 pc) stellar cluster 
as mapped by Thatte et al. (1997). 



4.4. The rotation curve 

The molecular gas velocity field inside the spiral arms is complex, making it hard to 
derive the rotation curve from ^^CO alone. We now discuss our best estimate of the CO 
rotation curve and compare our results with the stellar motions. We show that the velocity 
fields of gas and stars are similar, refiecting the infiuence of the bar on the motions at 
3" < r < 15". We think the difference in the rotation curves of the gas and stars could be 
due to the complex dynamics induced by the stellar bulge. 



4-.4-1- Rotation curve from gas motions 

The ^^CO velocity field in the center of NGC 1068 deviates strongly from circular 
rotation at r < 13", so a rotation curve derived from a rotating disk model would probably 
be wrong. The velocity field and pv diagrams show however, that most of the ^^CO emission 
is indeed from a centro-symmetric structure. We therefore used the first-moment velocity 
maps of the ^^CO (1-0) and ^^CO (2-1) emission to determine the velocity extrema in 
circular annuli around the center. We chose the ^^CO (2-1) velocity map because the 
high- velocity-dispersion knots on either side of the nucleus are partly resolved (Figs.[l| and 

The mean velocities in these maps are lower limits to the true rotation speeds since we 
did not correct for inclination. Our mean curve, corrected for an inclination i = 40°, agrees 
well with published rotation curves within the uncertainties (Fig. |5p. We compared our 
rotation curve with that derived from Ha (Dehnen et al. 1997; errors 10 to 30kms~^ for 
r < 16") and with the curve derived from the BIMA ^^CO (1-0) data (Heifer & Blitz 1995; 
errors 30 to 70kms-^ for r < 16" and 10 to 20kms-^ for 16" < r < 30"). 

From r = 2.3" moving inward, our ^^CO (2-1) rotation curve first increases by 50kms~^ 
and then decreases again at r = 1". This range of radii contains the nuclear ring and the 
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two barely-resolved CO knots east and west of the nucleus. The pv diagrams along p. a. 
110° are similar in ^^CO (1-0) and ^^CO (2-1) . We derived the rotation curve at r < 1" for 
two cases: (1) rotation velocities go to Okms^^ at the center, and (2) Keplerian velocities 
around a point mass interior to r = 0.4" (section |0| ). 



4-4-2- The stellar velocity field 

Garcia-Lorenzo et al. (1997) mapped the stellar velocities in the inner 24" x 20" with 
the optical Ca II triplet at A 8498, 8575, and 8695A. The velocity field of the stars matches 
that of the ^^CO (1-0) emission to within 50kms~^ and 1". This indicates that at r > 3" 
(the resolution of the Garcia-Lorenzo et al. data), the stars and gas have similar kinematics. 
The velocity dispersion observed in the Ca II lines agrees well with the velocity dispersion 
derived in the NIR (Dressier 1984; Terlevich et al. 1990; Oliva et al. 1995), suggesting that 
the velocity field measured in the Ca II line is indeed representative of the overall nuclear 
region (inner 15" to 20") and does not merely trace the surface of an otherwise obscured 
stellar system. 

For both stars and gas, at radii < 10", the isovelocity lines in the north bend to the 
east, and the isovelocity lines in the south bend to the west. Such a bending is as expected 
from our model calculations (section [7.2|) for orbits induced by a bar. It also agrees with the 
expected behavior of gas and stars in the potential of the NIR bar (e.g., Athanassoula 1992; 
Maciejewski & Sparke 1999). To confirm this finding observationally, more data are needed 
on the Ca II velocity field south (~ 8") of the nucleus, since the bending in the isovelocity 
lines occurs at a higher velocity relative to systemic than in the north. Higher-sensitivity 
molecular line data south of the nucleus would also be useful. 



4.4-3. Rotation curve in the inner Q" 

In the inner region the kinematics of stars and gas may be decoupled. The stellar 
kinematic data (Garcia-Lorenzo et al. 1997) indicate that the stellar rotation curve lies well 
below the rotation curve derived from the ^^CO data. This probably reflects the fact that 
in the inner 6" the velocity dispersion of the stars is higher relative to that of the molecular 
gas. In this region, the gravitational force of the gas is negligible relative to that of the 
stars, so we estimate the effects of different stellar potentials as follows. (1) We start with 
the simplest model, an isotropic, spherical stellar cluster. (2) We then generalize our model 
to central oblate stellar systems. (3) Finally, we consider the effect of a stellar bar on stellar 
orbits in the center of NGC 1068. 

Isotropic, spherical stellar cluster: The mass of an isotropic, spherical, non-rotating 
stellar cluster can be estimated from the stellar velocity dispersion and the Jeans equation, 
derived from the coUisionless Boltzmann equation (see Binney & Tremaine 1987). Thatte 
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et al. (1997) assumed a velocity dispersion of 140kms~^ (Dressier 1984) to estimate a mass 
of 6.5 X 10^ Mq within the inner 2" diameter. This is nearly twice our value of 3.7 x 10^ M© 
for the same region derived from the molecular gas. The discrepancy could be due to (1) 
the stellar flux coming from supergiants with deeper Ca II absorption than in the template 
giant stars, yielding a lower apparent velocity dispersion, or (2) the velocity dispersion of 
the underlying bulge component, which emits up to 45% of the light in the inner 5" (Table 
1 of Thatte et al. 1997). Comparison of /-band data in the inner 15", which show constant 
velocity dispersion (Dressier 1984), with the nuclear H and K band data indicates that the 
velocity spread of the nuclear cluster cannot exceed that of the outer zone. All this suggests 
the mass derived by Thatte et al. (1997) is an upper limit to the enclosed nuclear mass. 

Central oblate system: Observations and theory suggest bars turn into bulges at small 
radii (Norman et al. 1996; Hasan et al. 1993; Combes et al. 1990). For these systems 
the Jeans equation separates into two equations giving the vertical and radial form of the 
gravitational potential (see review by Merritt 1999). Binney & Tremaine (1987, their Fig. 
4-5) show that if self-similar isodensity ellipsoids of stellar spheroidal systems become more 
spherical at smaller radii then v/a falls below unity. This is indeed indicated by i^'-band 
images for r < 5" (e.g.. Fig. 1 in Thatte et al. 1997). Hence if the observed stellar velocity 
dispersion is constant, the circular velocity must decrease at smaller radii. This is consistent 
with the observed drop in the stellar rotation curve toward the center of NGC 1068. 

Influence of a bar potential: NGC 1068 has a stellar bar in the inner 20" seen in the 
NIR (Scoville et al. 1988; Thronson et al. 1989). Bars have prolate symmetry and are 
probably triaxial. Such a system can have an anisotropic velocity dispersion and its ellipsoid 
need not be coaxial with the bar (Fillmore 1986). It is thus not easy to predict the velocity 
dispersion, since in a bar potential, many orbit families are extended in the 2;-direction and 
there can be radial and vertical resonances that support specific orbit families, affecting 
the stellar velocity dispersion (e.g. Pfenniger 1984; Combes et al. 1990; Olle & Pfenniger 
1998). Unlike the stars in the center of NGC 1068, the molecular gas has a low velocity 
dispersion (section |4.1|) , indicating it is decoupled from the stars. This explains why the gas 
rotation curve — which we used for our models — lies above the stellar rotation curve by 
Garcia-Lorenzo et al. (1997) at small radii. Hence, for circular orbits, the molecular gas is 
a good tracer of the potential and its rotation curve is governed by the enclosed dynamical 
mass. 



4.5. The bars and resonances 

The interaction of the molecular gas with the stellar gravitational potential leads to 
resonances that in turn determine the distribution and kinematics of the gas. We now 
discuss the evidence for bars and the locations of dynamical resonances in NGC 1068. 
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4-5.1. The Bars 

Although NGC 1068 is not normally considered a barred spiral, Scoville et al. (1988) 
nevertheless detected at K band a 32" long bar at p. a. ~ 48°. The ends of the NIR stellar 
bar join the spiral arms in the molecular gas (e.g., Tacconi et al. 1997). In most bar 
theories, there are no massive pile-ups of gas at the corotation radius near the ends of the 
bar. These accumulations are only expected at the inner and outer Lindblad resonances. 
The fact that the arms are at the bar corotation radius is thus puzzling and could indicate 
the dynamical effect of a second, larger bar in NGC 1068. (Here, we assume that the 
corotation lies about at the end of the bar as found for late type galaxies, see Combes 1997 
for a review.) 

Optical images of NGC 1068 show a bright inner disk at radii < 100", with spiral 
arms and dust lanes extended north-south (see sketch in Fig. ^). This oval structure is 
180" (~ 13kpc) long, at p. a. 5°, with an eccentricity of 0.8 (Bland-Hawthorn et al. 1997, 
from BRI images by R.B. Tully). At radii > 15", however, the kinematic major axis runs 
almost exactly east-west. Hence the eccentricity of the oval structure on the sky is not due 
to the galaxy's inclination; instead, this structure is a large-scale bar along the minor axis 
of NGC 1068. In the Palomar Sky Survey images of NGC 1068, the outer, weak emission 
(radii > 100") is extended east-west, and the eccentricity of its lowest contours corresponds 
to the 40° inclination derived from HI velocities (Brinks et al. 1997). Deprojection of the 
optical image then yields a bar length of 240" (17.3 kpc). As shown in section |4.5.2| , this 
bar has its inner Lindblad resonance (ILR) at r = 18" (1.3 kpc), near the molecular spiral 
arms. The outer bar-like structure in NGC 1068 could thus be the cause of these spiral 
arms. Further support for this idea comes from the deprojected shape of the molecular 
spiral arms, which is an almost perfect ring — as expected for an ILR. It suggests the NIR 
bar is a second bar inside the ILR, similar to nuclear bars in other galaxies (e.g. Friedli et 
al. 1996). 



4.5.2. Location of the Resonances 

We estimated the radii of the dynamical resonances of the inner and outer bars from 
the rotation curve and bar lengths, using the curves of Q and Q ± k/2, to see if they 
matched the sites of the gas (Fig. |^ . (We follow the usual nomenclature, with Q = angular 
velocity, and k = epicyclic frequency, where = r x (dQ/dr)'^ + AQ"^ = 2Q x [Q + (dv/dr) 
— see, e.g. the review by Sellwood & Wilkinson 1993). 

We approximated dv/dr with the gradient Aw/Ar, which means the derived curves of 
K and (fi ± k/2) can have errors due to the undersampling of the rotation curve. The mean 
error in the angular velocity at radii > 5" is ± 12kms^^/kpc and larger for radii < 5", 
where the velocity field is more complicated. The local maximum in the {VL — /t/2) curve 
at r ~ 17" corresponds to the maximum in the rotation curve. The change in the velocity 
gradient at this point results in a change in k, which yields the maximum in {Vl — k,/2). 
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The change at this radius also occurs in the angular velocity curves of Heifer & Blitz (1995; 
their Fig. 18) and the stellar curves from Telesco & Decher (1988; their Fig. 6). 

The outer or primary bar has a deprojected radius of 120" (8.7 kpc). At this radius, 
the HI rotation velocity is 165kms~^ (Brinks et al. 1997), so the pattern speed of the 
primary bar is flp^ ~ 20kms~^/kpc with an uncertainty of ±10kms~^/kpc because the 
end of the bar does not exactly coincide with the corotation radius. As shown in Fig. |^ 
this pattern velocity Qp^ yields an outer ILR (oILR) at a deprojected radius of 1.6kpc (17" 
on the sky) and an inner ILR (ilLR) at 1.0 kpc (11" on the sky). This is identical with the 
location of the gas spiral arms at a deprojected distance of 1.3 kpc (~ 14" on the sky). It is 
thus likely that the gas spiral arms are caused by the outer bar driving molecular gas into 
the resonance zone, where the gas becomes trapped. 

If an inner bar exists, it should have its corotation radius at the ILR of the outer bar. 
This means the inner, NIR bar should have a deprojected radius of ~ 1.4 kpc, as observed. 
This radius yields an inner-bar pattern speed of Qp^ ~ 140kms~^/kpc, which allows an ILR 
of the inner bar at r < 6" — in the zone where we observe the nuclear gas ring. At these 
small radii, however, we know the rotation curve and {Q — k/2) only poorly. In contrast to 
our ideas. Heifer & Blitz (1995) estimated a deprojected NIR bar length of only 15" and 
thus derived an pattern speed of of 160kms~^. In their interpretation, the molecular spiral 
arms lie between the outer Lindblad resonance and the corotation radius. 



5. MOLECULAR GAS MASS AND DYNAMICAL MASS 

We now discuss the molecular gas and dynamical masses for some components of 
the ^^CO line emission. These values are used to estimate the thickness of the gas disk 
(see section ^ and the torques acting on the disk, should it be warped (see section |^ and 
Appendix C in SET99b). Comparison of the molecular gas mass with the dynamical mass 
of the nuclear region shows that the molecular gas is a probe of the nuclear gravitational 
potential in NGC 1068, but not a dominant component of that potential. 

Since the interferometer maps may have missed up to one half the total line emission, 
we derived masses only for the compact components (section ^4.2|) . The measured fluxes 
5*00^^ [Jykms~^] were converted into beam-diluted, integrated brightness temperatures, 
in Kkms~^, and then multiplied by the Milky Way conversion factor of 2 x 10^° 

j^'^™ from Strong et al. (1989) to estimate H2 column densities. We then multiplied by 
the beam-broadened source areas, in cm^, to obtain H2 masses (Table The calibration 
errors are only 15% (in the ^^CO (1-0) line), so the main uncertainty in the masses comes 
from the -^^-conversion factor, which may be wrong by at least 50% in the compact 
components, and even more in the extended components. In these mass estimates, since the 
integrated [S^qI^'qI line ratio is close to unity, we assumed the two CO lines were equivalent, 
provided maps were convolved to the same resolution and intensities were converted to 
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Kkms ^. 

We estimate the dynamical mass, in Mq, from Mdyn = 232 v'^{r) r, where v{r) is the 
inchnation-corrected rotation curve, in kms~^, and r is in pc. The uncertainties in this 
mass correspond to the uncertainties in the rotation curve (section ^.4\ ). Our value for the 
mass within r = 2" (for circular orbits) agrees with the upper limit from Thatte et al. 
(1997). Most of the gas in the inner 50" of NGC 1068 is in the molecular spiral arms, and 
its mass is comparable with that in the molecular ring of the Seyfert galaxy NGC 3227 
(SET99b). The molecular gas mass is only 5% of the dynamical mass in the inner 50". 
Even if one corrects for CO emission resolved out by the interferometer, the gas mass is less 
than 10% of the dynamical mass. 



6. DISK HEIGHT FROM THE VELOCITY DISPERSION 

The velocity dispersion observed in the molecular line emission (see Fig. and d) 
is due to (1) superposition of gas at different locations on the line of sight with different 
velocities in the same beam and (2) the intrinsic turbulent velocity dispersion of single 
molecular clouds. It is this turbulence that supports the disk and determines its height. 



6.1. The velocity dispersion 

The velocity dispersion was measured in the second-order moment maps of the ^^CO 
lines. In the regions of enhanced flux density on the spiral arms the velocity dispersion rises 
to 20kms~^ whereas the mean value in the spiral arms is ~ 16kms~^. This mean value is 
similar to the HI velocity dispersion of lOkms"^ for radii > 3 kpc (Brinks et al. 1997). A 
higher HI dispersion of 30 to 50kms~^ is observed at smaller distances from the nucleus 
(Brinks et al. 1997). The discrepancy with the velocity dispersion measured in ^^CO may 
be due to the larger HI beam (8") at radii < 20", or to the HI disk being thicker than the 
molecular gas disk, as in the outer regions of our Galaxy (Malhotra 1994, 1995). For the 
bar region we also find a low dispersion of ~ 16kms~^ in the molecular gas. 

In the nucleus, the mean velocity dispersion measured in ^^CO (2-1) is 40kms~^, twice 
as high as in the spiral arms. This may be due to inclination or beam smearing rather than 
intrinsic turbulence in the gas. A better estimate of the velocity dispersion of the nuclear 
gas may be the 25 km s~^ spread in the systemic velocity component in the pv diagram 
at p. a. 110°. This is probably the value that sets the disk thickness. The high apparent 
velocity dispersions in the two nuclear emission knots are due to the superposition of at 
least two distinct components (section [4.2|) . 
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6.2. The disk thickness 

To estimate the disk height from the velocity dispersion we must correct for beam 
smearing. We assume the observed dispersion dobs is the quadratic sum of the intrinsic 
dispersion, o"reai, and a contribution from the rotation curve (Xrot- We adopt an apparent 
velocity dispersion of (Xobs.sp ~ 16kms~^ for the spiral arms and aobs.nuc ~ 25kms~^ for 
the nuclear region (see above). At the spiral arms, the rotation curve is almost fiat, so in 
this zone, the CO line widths are not affected by beam smearing. In the nuclear region, 
for 0.5" < r < 1.0", the gradient in the rotation curve reaches 24kms~^/arcsec, which in a 
0.7" beam contributes cTrot,nuc ~ lOkms"^ . This only slightly increases the total apparent 
velocity dispersion o"obs,nuc- After deconvolution, the intrinsic velocity dispersion in the 
nuclear region is CTreai.nuc ~ 23kms~^. 

There are several ways to relate the velocity dispersion a and disk height h. Assuming 
hydrostatic equilibrium they can be written in their basic form as h/r ~ cr/uc, with Vc being 
the circular velocity at a radius r. Quillen et al. (1992) and Combes & Becquaert (1997) 
derived relations for gas disks in the potentials of elliptical galaxies or galactic bulges. 
Downes & Solomon (1998) deduce a relation for disks with flat rotation curves following 
the Mestel (1963) approximation. Both approaches are relevant to the center of NGC 1068, 
where the potential changes from disk-like to bulge-like. 

For the spiral arms, with a rotation velocity f(18") ~ 200 km s^^ and a dispersion 
Creai.sp ~ 16kms~^, the equation 3.1 of Quillen et al. (1992) yields a disk height of 100 pc. 
For the nuclear region, with v{l") ~ 150kms~^ and o"reai,nuc ~ 23kms~"'^, we get the 
following estimates. With eq. (3.1) of Quillen et al. (1992) we get a gas disk thickness of 
10 pc. With the (third) equation of Combes & Becquaert (1997), we find a disk height of 
9pc. The application of the Mestel approximation as given by Downes & Solomon (1998, 
their eq. 3) with areai,nuc ~ 23kms-\ w(3") ~ 170kms~^ and M(H2) ~ 5.3 x 10^ Mq yields 
an estimate of 7pc. All three methods thus indicate the nuclear gas disk is thin {d ~ 10 pc). 



7. DYNAMICAL MODELING 



As noted in section [4.1| , the velocity field of the central 2.5 kpc of NGC 1068 deviates 
from that of a simple inclined rotating disk. An obvious deviation occurs at r ~ 10", where 
the isovelocity contour at the systemic velocity, which should be a straight line along the 
kinematic minor axis, is in fact inclined by 45° to the kinematic minor axis of the outer 
disk. Another difference to a simple disk is the complex structure in the velocity field and 
pv diagrams of the inner few arcseconds. Our kinematic models can fit these structures by 
either planar elliptical orbits due to gas moving in bar potentials, or circular orbits tilted 
out of the galactic plane in a warp. We first describe the best fits in pure bar models and 
pure warp models, and present a combined solution as the most plausible explanation. We 
then relate this best-fit, mixed model to observations at other wavelengths. 
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7.1. Outline of our model program, SDRings 

Details of our SDRings program are in SET99b (appendix B), so we give here only a 
brief summary. We assume that since gas is dissipative, it persists only on orbits that are 
not self-intersecting, not crossing, and not strongly cusped. Its orbits can thus be of only 
two basic types: (1) planar elliptic orbits and (2) tilted circular orbits. The first type are 
the Xi and X2 orbits along and perpendicular to bars in the galactic plane. The second type 
are orbits that can leave the galactic plane, thereby warping the gas disk. 

The model divides the disk into many individual circular or elliptical orbits of molecular 
gas. For the modeling the inclination and position angle of the disk, and the shape of the 
rotation curve were held fixed. Each fit was started at large radii and then extended to the 
center. For each case we tried several setups that all converged to similar best solutions 
with mean deviations from the data of < lOkms"^ and 0.1" for each velocity and radius in 
the pf-diagrams and 10° in the position angle of the mapped structures. 

For the elliptical orbit models, we follow Telesco & Decher (1988), who explained the 
gas motions between the two ILRs in NGC 1068 with a smoothly varying position angle of 
the ellipses centered on the nucleus. In the bar models, we fitted the orbital eccentricity 
e(r) = b{r)/a{r) (where a and b are the major and minor axes) and the position angle 
P Ae\\[psc,{r) by curves varying smoothly with radius, with the orbits constrained to be 
non-overlapping. The orbits lie in a plane and resemble the velocity and density patterns in 
the bar models of Athanassoula (1992). 

For the out-of-plane circular orbits, we follow the tilted ring model outlined by Rogstad 
et al. (1974) with a warp simulated by concentric tilted rings. The strength of the tilt and 
the precession of the position angles are proportional to the acting torque. The torque can 
result from gas moving in a non-axisymmetric potential, or from interaction of the disk gas 
with the radiation pressure of a nuclear source, with the gas pressure in the ionization cone 
of an AGN, or with some GMCs out of the disk (see SET99b, Appendix C). 

To model the warp, we followed the method of Tubbs (1980; see also Quillen et al. 
1992), in which the warp has a smoothly varying tilt a;(r) of each orbit relative to the disk 
plane and a smoothly varying precession angle a{r) (Fig. ^). A torque acting on an orbit 
with a circular velocity Vc{r) induces a precession rate da/dt ~ ^Vc/r. After a time At one 
obtains a{r) = ^fiAt + ao, where ao is a constant, Q=Vc/r, and ^ is given by the acting 
torque. Our models had constant ^At and uniformly distributed gas. 

The model fits were adjusted to match the observed pv diagrams, the velocity moment 
map, and the intensity map. Comparison of the models with the data was done in maps at 
a resolution of 0.4" which is super-resolved relative to the synthesized beam . This approach 
allowed us — especially in pf-diagrams — to concentrate on the brightest compact features 
in the data with a 0.7" beam. For both approaches we tried two rotation curves, one 
interpolated to Okms~^ at r = 0" and one in which the observed rotation curve in the inner 
few parsecs was replaced with a rotation curve of a fixed enclosed mass. The best results 
were for a warp model with an enclosed mass of ~ 10^ Mq. 
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7.2. The pure bar approach 

We showed above that two bars are clearly detected in NGC 1068 and that the spiral 
arms are at the ILR of the outer bar (see section |4.5| ). In the inner 4" the potential is 
dominated by the inner 1 kpc NIR bar. To be stable, an inner bar must be weaker than an 
outer one (e.g., Friedli & Martinet 1993). The influence of the NIR-bar on the gas motions 
is thus weaker than that of the outer one (Macjiewski & Sparke 1999). Furthermore, if the 
inner bar controlled the gas motions, the large velocity dispersion observed at a radius of 1" 
would be hard to explain, since strong streaming motions are not expected for inner bars. 
There is also no hint in the NIR data (Thatte et al. 1997) for an additional, strong, third 
bar that could influence the structure and kinematics of the circumnuclear molecular gas. 

For our pure bar model. Fig. ^ shows the model curves for the position angle and 



ellipticities of the fltted ellipses, with their errors (see section |7.1| ). Figure [T0| compares the 



model pv diagrams with the data, and Fig. 11 shows the model intensity distribution and 



velocity fleld. We now discuss the different regions of the model. 

Spiral arms and bar region: Our best bar model gives a good flt to the data cube 



in the spiral arms and in the transition to the NIR bar (Fig. [Til) . In particular, the 
model reproduces the shifts in the alignment of the kinematic minor axis. Inspection of the 
velocity fleld on the northern part of this axis reveals the behavior of the molecular gas at 
an ILR and near the inner bar. At an ILR (the spiral arm region in NGC 1068), the velocity 
is blueshifted relative to systemic, and becomes redshifted near the NIR bar. That is, at 
the 30" diameter spiral arm region, the velocity on the minor axes in both ^^CO (1-0) and 
^^CO (2-1) is blueshifted relative to systemic because of the ILR of the outer large-scale 
bar. This velocity change along the minor kinematic axis occurs close to the spiral arms, at 
the transition from the Xi orbits (those running along the outer bar, outside the oILR) to 
the X2 orbits (those normal to the outer bar, between the oILR and the ilLR). This outer 
velocity change is due to the X2 orbit velocities being high where the xi orbit velocities are 
low. 

At radii < 7", the velocity is redshifted because of the NIR bar. At the ilLR, there is a 
transition from X2 orbits of the outer bar to xi orbits of the inner, NIR bar at p. a. 48° 
(Thronson et al. 1989). This is the second, inner, velocity change along the kinematic 
minor axis. In NGC 1068 this second shift is seen in the ^^CO (1-0) and ^^CO(l-O) lines 
(see section ^ and Tacconi et al. 1997; Heifer & Blitz 1995). This scenario is well matched 
by the orbit crowding in our model. 

The nuclear gas ring: Near the nuclear ring, the bar model fails to reproduce the 
observed velocity increase moving in to 1", the subsequent abrupt drop to the systemic 
velocity, and the observed ridge between the two CO knots (section Fig.^. This is a 
problem for the pure bar model, and one may ask if it could be rescued with a third bar. 
Fitting of the gas motions with the bar approach suggests the presence of highly elliptical 
east-west orbits and therefore provides evidence for the the presence of a third bar with a 
major axis of 1" length. High-resolution NIR observations, however, show no such nuclear 
bar in NGC 1068 (Wittkowski et al. 1998; Weinberger et al. 1999; Thatte et al. 1997). The 
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stellar continuum in the nucleus has neither the 1" size nor the shape expected for a third 
bar (see Fig.2 in Thatte et al. 1997). 



7.3. The pure warp approach 

Another model for non-circular motions is a warp in the gas disk. The warp model 
yields a solution that is fully consistent with the derived gas rotation curve. In particular, it 
can explain the high velocities between 1" and 2" in the nuclear ring, seen in the pv diagram 
at p. a. 110°. (As in the bar model, for the inner few parsecs we adopted a Keplerian 
rotation curve for an unresolved point mass — cf. section |7.1| .) 



The precession of the warp is given by the torque ^, the time interval At, the initial 
precession angle cto, and the angular velocity Q{r) of the rotation curve. For the best 
model, we get ao = 10° ± 15° and ^At = (7.0 ± 0.5) x 10^ yr. In this model, it is the fit to 
the intensity distribution that is more sensitive to changes of ^At, not the fits to the pv 
diagrams. The precession time is much longer than the rotation time at r = 1". At this 
radius, with our preferred value of ^At, the disk must make 25 rotations in order to precess 
by 360°. 



The inclination curve uj^r), which is the tilt of the warp, is shown in Fig.[T2| for the 
best fit. The errors range from ± 5° to ± 30°. The small errors at r = 1" are due to the 
abrupt drop in the systemic velocity at this radius, which in turn depends strongly on the 
inclination of individual rings. 

The spiral arms and bar region: A warp can mimic spiral patterns in the plane of the 
sky (Steinman-Cameron et al. 1992). What about the molecular spiral arms in the center 
of NGC 1068? Their kinematic axis changes direction, suggesting these arms depart from 
circular rotation and may also be a warp effect. We found that to fit the spiral arms, the 
model disk must change its inclination by 20° between r = 20" to 17", then stay at this 
new inclination down to r ~ 2.6", and then flip back to the old inclination at r = 2.2". One 
anomaly is the kinematic minor axis, which in the observed data, crosses the zero-velocity 
contour twice, at the spiral arms. The warp model cannot reproduce this effect. Otherwise, 
the warp model fits the velocity field inside r = 17" just as well as the bar model does. 

For the nuclear ring, Fig. ^ compares the observed and model pv diagrams, and 
Fig. |l^ shows the intensity map and the velocity field of the warp model. For the symmetric 
pv diagram at p. a. 110° with most of the line flux, the model deviates from the data by only 
±5kms~^. In the nuclear ring, the model can explain everything in the pv diagram and the 



spectra (Fig. JSj): the large velocity dispersion, the associated increase in velocity at 1", 
the emission at the systemic velocity at r < 1", the decrease of velocity toward the center, 
and the ridge between the two high dispersion points. On the minor kinematic axis, the fit 
is not as good, with deviations of ±20kms~^. This may be due to the low line flux and its 
contamination by sidelobes that are not completely removed by the CLEAN algorithm. 
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For radii < 2.2" the disk starts warping, down to r ~ 1.5" (110 pc) where it becomes 
almost edge-on (i' = i + ~ 70°). It stays at this inchnation for about 35 pc (0.5") toward 
smaller radii, and then warps within r < lOpc from edge-on to almost face-on. Due to 
resolution limits, the orbit orientations at r < 30 pc are not well defined. The outer warped 
molecular gas disk (10 pc to 100 pc) probably joins continuously the inner accretion disk 
(~ Ipc), which may also be warped (Gallimore et al. 1997; Greenhill & Gwinn 1997). The 
warped molecular gas disk may thus be not only the link between the host galaxy and the 
accretion disk, but also the fuel reservoir for the AGN (see section 8). 



7.4. The mixed model 

The bar and warp approaches both provide good fits to the data for specific regions 
of NGC 1068. From the observations of the two bars and the analysis of the dynamical 
resonances, it is obvious that the modeling of the spiral arms with elliptic orbits representing 
a bar potential is plausible and physically realistic. This is not necessarily true for the 
nuclear region, since the potential of the inner stellar (NIR) bar is weaker relative to the 
outer bar, and the warp approach provides a better fit to the individual components. 

We thus use a hybrid model to describe the kinematics of NGC 1068, combining the 
best fits from both approaches. In the warp model, the tilt of the nuclear disk goes in the 
opposite direction of that needed to model the spiral arms. This allows us to join, at the 
outer edge of the nuclear ring, the warp model for the inner zone with the bar model for 
the outer zone. This hybrid model fits very well the intensity and velocity distributions of 



all the symmetric features. Figure |T6| is a perspective view of the 3-dimensional model in 
the plane of the sky. In this mixed model, the spiral arms are at the ILR of the outer bar. 
The NIR bar is an inner bar with its Xi orbits causing the change in the position angle 
inside the spiral arms. The bar model of both the ILR of the outer bar (at the spiral arms) 
and the transition to the inner bar fits well the observed velocity field. The nuclear ring, 
however, cannot be reproduced very well as an ILR of the inner bar, in the same plane. 
The warp approach gives a better fit to the nuclear gas ring. 



8. DISCUSSION 

In section [1.5| we described the 17kpc-long outer bar and the 2.7kpc-long, NIR inner 
bar. Gallimore et al. (1999) recently resolved the HI absorption in the nuclear region of 
NGC 1068 in both space and velocity. The strongest HI is 4" southwest of the nucleus and 
absorbs continuum emission of the southern radio lobe (Fig. 4 of Gallimore et al. 1999). 
The HI absorption velocity varies from southeast to northwest. Gallimore et al. (1999) 
modeled this motion with a disk inclined at the position angle of the stellar bar, and 
suggested the motion may be due to gas in the potential of that bar. The velocity fields of 
our ^^CO (1-0) data and our bar model confirm this conclusion. Comparison of systemic 
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velocity measurements (Heifer & Blitz 1995; Brinks et al. 1997; Greenhill & Gwinn 1997; 
this paper) suggests the value for NGC 1068 is ~ 1120 km s~^ rather than the 1150 km s~^ 
derived by Brinks et al. (1997) from their HI data with a 8" beam. This implies that the 
HI absorption at 1110 and 1131 kms~^, seen in projection near the nucleus in the maps of 
Gallimore et al. (1999), is actually physically associated with the nucleus. 

Friedli & Martinet (1993) made N-body simulations of double-barred galaxies and 
showed the dynamics inside the outer bar's ILR can be decoupled, so that an inner bar can 
form. This is consistent with our bar model for the ^^CO data. The inner and outer bars 
are usually assumed to have similar ILRs (Friedli & Martinet 1993; Maciejewski &; Sparke 
1999). So far, however, simulations of inner-bar ILRs have not been done with high enough 
numerical resolution to yield the details of such a region. 

For the outer bar to survive, the inner bar's potential must be much weaker, implying 
the inner bar has a less pronounced ILR than the outer bar (Friedli & Martinet 1993). The 
molecular gas ring at a radius of 75 pc from the nucleus of NGC 1068 is close to the central 
stellar cluster (Thatte et al. 1997). The gravity of the central mass concentration competes 
with that of the inner, NIR bar, further weakening the inner bar's influence on the gas in 
the ILR region. Near the inner bar's ILR, a vertical resonance may evolve (e.g. Combes et 
al. 1990). Since only non-crossing orbits can retain the gas, such a vertical resonance would 
warp the inner molecular disk. Once out of the disk, the gas would be further supported by 
the potential of the central stellar cluster, competing with that of the bar. 

Advantages of the warp model for the nuclear region: In summary, the reasons favoring 
a warp model are: (1) The best warp model reproduces all the observed symmetric 
structures in the CO pv diagram along p. a. 110°. (2) The warp forms an opening for the 
ionization cone (Fig. [T^). (3) The warp provides enough dust between the observer and 
the AGN to explain why the NIR peak is offset from the optical one (see below). (4) In the 
inner 3", the polarization vectors are aligned as expected for the edge-on geometry of the 
molecular gas in the warp model (see below; an edge-on geometry for the central few arcsec 
was also proposed by Baker & Scoville (1998) and Baker (1999) to explain their CO data). 

Cause of the warp: To move gas out of the plane, a torque is needed (see appendix 
C in SET99b). Torques can be induced by a non-spherical galactic potential (just as 
the halo acts on the HI disk), by radiation pressure from a radio jet (similar to central 
radiation sources causing warps in accretion disks; Pringle 1996; 1997), by gas pressure in 
an ionization cone (Quillen & Bower 1999; SET99a,b) or by the transient gravitational force 
of a giant molecular cloud temporarily dislocated above the plane of the disk. Arnaboldi & 
Sparke (1994) calculated the torque of a slightly oblate galactic potential outside the core. 
We assumed a similar potential, and used their eq. (B7) to estimate the torque. Our results 
(Table ^) show the most likely causes of a warp are torques induced by gas pressure, or 
by dislocated GMCs. The torque induced by the gas pressure of the ionization cone, for 
example, could persist and maintain the warp for ~ 10^ yr. A warp caused by a CMC (e.g., 
the western knot), deviating from the general rotation, would only perturb the galactic 
potential for as long as the CMC is above or below the molecular disk. 



- 19 - 



Optical dust lanes: The dust lanes in the center of NGC 1068 are obvious in the HST 
image with the F550M filter (Fig. 3 of Catchpole & Boksenberg 1997). To enhance the 
contrast, the authors subtracted a model of the nuclear starlight from the image. Their 
model was azimuthally symmetric, following two power laws. The model-corrected HST 



image nicely shows enhanced extinction running east-west. Figure compares this 
corrected image with the intensity distributions from the bar and warp models. This 
comparison favors the warp model, which predicts absorbing material at about the right 
locations. In a planar bar model, the extinction in front of the nuclear region can be 
explained only by postulating a very thick nuclear disk, which is inconsistent with disk 
thicknesses estimated from the velocity dispersion of the molecular gas. 

The near- and mid-IR nuclear polarization: The observed near- and mid-IR polarization 
vectors in the inner 3" are aligned east-west (Packham et al. 1997; Young et al. 1996). 
For Mie scattering in an edge-on disk, light should be polarized parallel to the disk. The 
reason is as follows: Bands of parallel polarization vectors are often observed in bipolar 
stellar outflows and can be reproduced by multiple scattering models (Bastien & Menard 
1988; Fischer, Henning, & Yorke 1996). These bands are often normal to the flow and occur 
especially in edge-on disks with high densities and particle sizes that are large relative to the 
wavelength of scattered light (Fischer et al. 1996). Besides strong forward and backward 
scattering, the Mie process also yields highly polarized emission, with the E vectors and 
the maximum polarized intensity both perpendicular to the direction of forward scattering. 
The light is strongly concentrated in two outflow cones, with the intense, back-scattered 
light illuminating the disk. This scattered light is then polarized parallel to the disk, over 
the total illuminated area (Fischer et al. 1996; Bastien & Menard 1988). Bands with 
parallel polarization thus occur most clearly in edge-on disks. 

From the near- and mid-IR polarization maps of NGC 1068, Young et al. (1996) and 
Packham et al. (1997) derive scattering cones at p. a. 30° and a "torus" size of 200 pc (3"), 
consistent with the orbit diameters in our warp model, for the orbits with large inclinations 
to the line of sight. Alternatively, the alignment of the NIR polarization vectors could be 
due to absorption of the central stellar continuum by elongated dust grains ahgned in a 
magnetic fleld by the Davis-Greenstein effect (Greenberg 1978). If so, then the magnetic 
fleld lines should lie in the plane of the absorbing gas disk, parallel to the polarization 
vectors. In this scenario, however, the mid-IR polarization vectors should be orthogonal to 
the NIR vectors, since the scattering efficiency is low at long wavelengths (unless the grains 
are very large), and the mid-IR flux is dust emission rather than scattered light. 



9. CONCLUSIONS 

1. Molecular gas close to the nucleus. — Our new interferometric data allow a 
high-sensitivity study of the molecular gas in the inner 1 arcmin of NGC 1068 at high 
angular (0.7" x 0.7") and spectral resolution (10kms~^). We made detailed maps of the 
200pc-diameter nuclear ring indicated in previous molecular-line maps (Tacconi et al. 
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1997). We now find indications of molecular gas emission as close as 13 pc from the nucleus 
with a large offset from the systemic velocity. This emission suggests a mass within the 
inner 25 pc of 10* Mq. This is consistent with the black hole mass of 1.7 x 10^ indicated 
by the H2O masers (Greenhill & Gwinn 1997), plus a significant contribution from the mass 
of the nuclear stellar cluster (Thatte et al. 1997). 

2. Bar and warp model. — Modeling of the molecular gas kinematics in the inner 50" 
suggests a hybrid model combining bar-induced motion in the spiral arms and NIR bar 
with a warp of the thin gas disk in the nuclear ring (inner 300 pc) . Due to the warping, the 
molecular gas leaves the plane of the disk at a radius of about 150 pc and has an inclination 
of ^ 90° to the line of sight at about 75 pc. 

3. Extinction of the nucleus. — The obscuration of the AGN by the warped molecular 
gas disk is consistent with the observations of the near- and mid-IR polarization and the 
observed extinction of the nucleus at optical wavelengths. It explains its identification as a 
Seyfert 2 galaxy with a hidden Seyfert 1 nucleus. 

4. Cause of the warp. — The probable cause for warping of the gas disk is the gas 

pressure in the ionization cone. This scenario is also proposed for the geometry in M 84 
(Quillen & Bower 1999). Further possibilities arc the influence of the transition between 
the disk/bar and bulge potential, or dislocated GMCs transiently disturbing the gas disk. 

5. Small tori are not needed. — As the results by Gallimorc ct al. (1997) show it is 
very likely that NGC 1068 has a <2 pc molecular torus. However, our findings imply that 
the nuclear molecular gas disks may be warped providing plenty of material to occult the 
nuclear regions. Therefore we do speculate that long-postulated (large-scale) molecular tori 
with radii of a few 10 pc to a few 100 pc are not needed in every case. The more likely 
explanation of the Seyfert 1/Seyfert 2 differences is that the molecular gas in the center 
of the host galaxy contributes significantly to the obscuration of the active nucleus and 
thereby influences the classification (Malkan et al. 1998). In our model, the absorbing 
molecular clouds need no longer be confined to the galaxy's plane, in agreement with the 
geometry already proposed by Cameron et al. (1993) for NGC 1068. This implies that the 
apparent nuclear properties are strongly influenced by the distribution of gas and dust in 
the center of the host galaxy. 

6. The gas spiral arms are at the ILR of a large-scale outer bar. — Analysis of the 
structure of NGC 1068 shows that the 3 kpc long, NIR stellar bar is an inner bar, since 
the galaxy also has an outer bar- like structure about 17 kpc long. This means that the 
spiral arms in the molecular gas (at a radius of 1.4 kpc) are at the ILR of the outer bar, in 
agreement with the theory of gas dynamics in barred potentials. Comparison of the velocity 
fields of the molecular gas and the stars with our model velocity field of a bar shows they 
are remarkably similar. 
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Fig. 1.^ Maps of (a) ^^CO (1-0) and (b) ^^CO (2-1) in the center of NGC 1068. Contours 
are at 10, 20, ... 100% of the peak, which is 1.24 Jy km/s/beam in ^^CO (1-0) and 3.43 Jy 
km/s/beam in ^^CO (2-1) . The central box in the ^^CO (1-0) map shows the area of the 
i^CO (2-1) map. 
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Fig. 2. — The distribution of the 110 GHz continuum in NGC 1068 with a sensitivity of la 
— 1 mJy/beam. The contours start at 3a in steps of 3a. The peak intensity is 23 mJy/beam. 
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Fig. 3. — Velocity field and velocity dispersion of the ^^CO lines in NGC 1068. The velocity 
field (a) and the velocity dispersion map (b) of the ^^CO (1-0) emission (beam 1.4") and 
the ^^CO (2-1) emission (beam 0.7") (c and d) are shown, respectively. The contours of the 
velocity fields are in steps of 20kms~^ (v=Okms~^ corresponds to the first central sohd fine 
next to a broken line). In the dispersion maps, contours are at 15, 30, 45 and 60kms~^ for 
^^CO (1-0) and at 15, 30, 45, 60 and 75kms-i for the ^^CO (2-1) . 
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Fig. 4. — Position-velocity diagrams of the nuclear ^^CO (2-1) line emission in NGC 1068 
along different position angles. To highlight this complex velocity behavior the data is shown 
at a nominal resolution of 0.4" (right). The contours are at 10, 30, ... 100% of the peak. For 
comparision the corresponding pv-diagrams at the achieved instrumental resolution of 0.7" 
are shown in color on the left panels. The velocity resolution in both cases is 10 kms~^. The 
(magenta) contour corresponds to 3cr, with lcr= 6.0 mJy/beam. A remarkable feature in the 
pf-diagram along p. a. 110° is the large velocity dispersion at ±1" from the center. For the 
first time in NGC 1068, we can trace molecular gas at ~ 0.18" (13 pc) from the nucleus. The 
high velocities in this region imply an enclosed mass of ~10^ Mq not correcting for inclination 
effects. This value is consistent with a black hole mass of 1.7 x 10'' M0, as estimated from 
nuclear H2O maser emission (Greenhill & Gwinn 1997), plus a contribution from a compact 
nuclear stellar cluster (Thatte et al. 1997). 
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Fig. 5. — Rotation curves for NGC 1068. The CO curves are the average values of the 
minimum and maximum velocity at any given radius (uncorrected for inclination). The other 
curves are taken from the quoted references. The increase in velocity of the ^^CO (2-1) curve 
at a radius of r ~ 1" can also be seen in the pv diagrams of the observed data. At radii 
between 10" and 16" the rotation curve agrees with fits of a rotating gas disk to the observed 
data. 
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Fig. 6. — Sketch of NGC 1068 showing the location of the bars and molecular spiral arms. 
With the exception of the outhne of the galaxy disk the figure is approximately to scale. 



-32- 




Q I , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 1 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 

2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 

Radius ["] 



Fig. 7. — Sites of inner Lindblad resonances (ILRs) and corotation resonances (CRs) in 
NGC 1068. The ILRp of the outer, primary bar is at the molecular spiral arms at r ~ 18". If 
an inner, secondary bar exists, its CR^ coincides with the ILRp of the outer bar. The outer, 
primary bar has a deprojected radius of about 120" translating into a CR at about 140" . 
(See text.) 



-33- 




Fig. 8. — Parameters in the SDRings warp model. The diagram shows only a single ring, 
tilted by an angle lo relative to the galaxy's rotation axis, which is inclined at angle i to the 
line of sight. Here a is the precession angle between the projection of the line of sight on 
the plane of the galaxy and the ring axis. 
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Fig. 9. — The curves of position angle PA and ellipticity e of the elliptical orbits in the bar 
model that fit the data. The curves of e (solid line) and PA (broken line) are shown as thick 
lines for the best fit, and as thin lines over the range of satisfactory fits (see text). 
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Fig. 10. — Position-velocity diagrams for the bar model in NGC 1068. The regions of 
high velocity dispersion that appear at low intensity (p. a. 10°, radius —1.8" and p. a. 150°, 
r = 1.3") are molecular cloud complexes that are not part of a symmetric velocity field (see 
section |4.2|) . The bar model explains all symmetric structures. Although the bar model can 
fit the component at lowest velocities, it cannot explain the high velocity dispersion of these 
components. The data at 0.4" resolution are shown in black contours of 10 to 100% in steps 
of 10%. The model, calculated for a starting radius of 27 pc, is shown in red contours at 0.2, 
2, 5, 20, 40, 60 and 80 %. 
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Fig. 11. — Intensity maps (left) and velocity fields (right) of the bar model for the ^^CO 
emission (a and b: ^^CO (1-0) , c and d: ^^CO (2-1) ). Velocity contours are from —150 
to 150 km s-^ in steps of 20kms-^ for ^^CO (1-0) and from -130 to 130 km s"^ in steps of 
20kms~^ for ^^CO (2-1) . Broken lines indicate negative velocities. 
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Fig. 12. — The u (angular velocity) curve of the warp model. The solid line shows the best 
fit to the data. Broken lines indicate the range over which the fits are satisfactory. 
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Fig. 13. — Position-velocity diagrams for the warp model in NGC 1068. The weak emission 
regions of high velocity dispersion (at p. a. 10°, r = —1.8" and p. a. 150°, r = 1.3") belong to 
molecular cloud complexes which are not part of a symmetric velocity field (see section |4^) . 
Our model only explains symmetric structures. The contrast-enhanced data are in black 
contours of 10 to 100% in steps of 10%. The model, calculated starting from r = 27 pc, is 
shown in red contours at 0.2, 3, 20, 50 and 80 %. 
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Fig. 14. — Intensity maps (left) and velocity fields (right) of the warp model for the ^^CO 
emission [top: ^^CO (1~0) , bottom: ^^CO (2-1) ). Velocity contours are from —150 to 
150 km s-^ in steps of 20kms-^ for ^^CO (1-0) and from -130 to 130 km s"^ in steps of 
20kms'"^ for ^^CO (2-1) . Broken lines indicate negative velocities. 
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Fig. 15. — Comparison of measured spectra (black) with those given by the warp model 
(red) and the bar model [green), in the central 2" x 2" in steps of 0.5". Intensity scales are 
arbitrary. With the exception of the eastern knot (knot E) the warp model represents an 
acceptable fit to the observed data that reflect the complex velocity field. The agreement 
of the bar model with the data is not as good as it cannot reproduce the multiplicity (two 
features) of the line, especially at positions west and south-west of the nucleus. 



Fig. 16. — Geometry of the NGC 1068 warp model. Brighter sections are closer to the 
observer, darker ones are farther away. The brightest structure is the [O III] ionization cone 
(Macchetto et al. 1994). The 3-dimensional geometry of the warp creates a natural cavity 
for the ionization cone that is consistent with its observed orientation. 
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Fig. 17.— The F550M filter HST image {top left) of tlie central 4" x 4" in NGC 1068 
divided by a model of the stellar light (see text and Catchpole & Boksenberg 1997). The 
dust lane is encircled by the solid gray line in {top right). The bottom right image shows 
that the CO intensity predicted by the warp model agrees with the dust distribution. The 
agreement is not as good for the bar model {bottom left), which does not predict the gas and 
dust concentration above and below the nucleus. 



Table 1: Properties of NGC 1068 



NGC 1068 

Right ascension (J2000) 02'^ 42"^ 40.798*^ 

Declination (J2000) -00° 00' 47.938" 

Classification (R)SA(rs)b 

Incfination 40° 

Position angle 278° 

AGN Type Sey 2 

Systemic velocity 1150 kms"^ 

Distance 14.4 Mpc 



Nuclear coordinates arc from Muxlow et al. (1996). Systemic velocity and Seyfert type are 
from NED (NASA/IPAC Extragalactic Database). The host galaxy classification is from de 
Vaucouleurs et al. (1991). Inclination and position angle are from the "Ringberg standards" 
(Bland-Hawthorn et al. 1997). 
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Table 2: Molecular gas masses of various compact components in NGC 1068 



component Sco M(B.2) Mdyn 

[Jy km/s] [10^ Mq] [IQS Mq] 

670 68 120 

560 57 

2140 5 

20 5 9 

12 0.3 

8 0.2 



total 

spiral arms 
Northern bar 
ring in ^^CO (2-1) 
knot E 
knot W 



Uncertainties in gas masses are 40% for the compact components and 50 % for the extended 
components. Uncertainties of dynamical masses are dominated by rotational velocity errors. 
The fluxes were derived for the total line width with exception of the two knots in the ring 
for which we used only velocities from —230 to —30 km s^^. The flux of knot E was measured 
in a rectangle with corners at —1.3"/ — 2.0" and 0.6"/ — 0.9" and knot in a rectangle with 
corners at 0.8"/ — 1.3" and 1.9"/ — 0.4". The dynamical mass was derived with the equation 
in the text, with our molecular gas rotation curve at radii of 18" and 2". 
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Table 3: Causes for the warp in NGC 1068 



Cause 


Property 


Value 




Estimated via 








SDRings" 


torque^ 


M 


1.9x10^^ Nm 




molecular gas mass 


m 


2.0x10^ Mq 




radius 


r 


2.2x10^8 m 




velocity 


v{r) 


140 km s"^ 




time of circulation 


$-1 


4.0x10^3 g 


potential'' 


torque^ 


M 


7.6x10^" Nm 




molecular gas mass 


m 


2.0x10^ M0 




volume 


Vo 


1.5x10^6 m^ 




total mass 


rrio = PoVo 


4.3x10^ 


GMC^ 


torque'^ 


m 


3.2 xlO^* Nm 




cloud mass 


^GMC 


2.0x10^ Mq 




force 


F 


9.5x1025 N 




lever arm 


I = r 


3.3x10^8 m 



see caption of next table 



NOTE TO EDITOR: The following two tables should be one; table caption is incomplete 
in processed file (ps-file) 
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Table 4: Causes for the warp in NGC 1068 continued 



gas pressure'' 


torque ^ 


M 


1.6x10^^ Nm 




particle density 


N 
V 


lO'' cm~^ 




temperature 


T 


1.0 X 10^ K 




gas pressure 


P 


1.4 X 10-^ Nm-2 




area 


A 


5.5 X 10^^ m^ 




lever arm 


I 


2.1 X lO^^m 


radiation pressure^ 


torque^ 


m 


4.8 X lO^'^Nm 




spectral index 


a 


~ -1.0 




constant 


h 


1.0 X lO-i^Wm-2 




luminosity 


L 


6.6 X 10^2 W 




intensity 


I 


1.2 X lO-^Wm-2 




radiation pressure 


P 


3.9 X 10-1^ Nm-2 




area 


A 


5.5 X lO^^m^ 




lever arm 


I 


2.1 X 10^8 m 



Physical reasons for the warp of circum-nuclear gas disks are summarized in the conclucions 
(section!^) and given in detail in the appendix C in SET99b. In this table first order estimates 
are based on: 

" from SET99b, eq. C3. ; ^ from SET99b, eq. C2. ; " from SET99b, eqs. CI and C4. ; 

from SET99b, eqs. CI and C5. ; " from SET99b, eqs. CI and C6. 

^ from the parameters of SDRings for an inner disk, with r = 108 pc 

^ from the mass of warped nuclear gas disk and a dynamical mass within r < 108 pc. 

^ from mass of knot W ai r = 108 pc and its distance relative to the gas disk. 

^ N/V and T of ionization cone from Capetti et al. (1997a) and Gallimore et al. (1996a,b) 

for jet components C and NE. A = 0.5" x 0.25" at r = 1", from radio jet map. 

^ spectral index, constant b, and area (0.5" x 0.25") of jet at r = 72 pc from Gallimore et al. 

(1996b). 



